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The band structure of silicon is calculated at the V, X, and L wave vectors using diffusion quantum Monte Carlo 
methods. Excited states are formed by promoting an electron from the valence band into the conduction band. 
We obtain good agreement with experiment for states around the gap region and demonstrate that the method 
works equally well for direct and indirect excitations, and that one can calculate many excited states at each wave 
vector. This work establishes the fixed-node DMC approach as an accurate method for calculating the energies 
of low lying excitations in solids. 
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I. INTRODUCTION 

Electronic excitations play a crucial role in the physics 
and chemistry of atoms, solids and molecules. Calculat- 
ing excitation energies in large systems is a challenge for 
theoretical techniques because an accurate description re- 
quires a realistic treatment of the electron correlations. 
The well-known Hartree-Fock (HF) method includes ex- 
change but not correlation effects, and therefore overes- 
timates band gaps and band widths by a large amount, 
while Kohn-Sham density functional calculations within 
the local density approximation (LDA) underestimate 
band gaps. Here we report a study of excitation energies 
in bulk silicon lisdpg the diffusion quantum Monte Carlo 
(DMC) methodtlli The DMC method is very promising 
for applications to condensed matter because (i) it ex- 
plicitly includes electron-electron correlation effects, and 
(ii) it scales reasonably well with system size, with the 
computational cost increasing as the cube of the number 
of electrons. 

It is now well established that DMC calculations can 
give an excellent description of electron correlations in 
the ground state. The range of problems that could be 
addressed using DMC would be greatly increased if one 
could also obtain accurate excitation energies. Further- 
more, the DMC method should be equally applicable to 
both strongly and weakly correlated systems. However, 
calculating excitation energies in condensed matter sys- 
tems is a formidable challenge to DMC techniques be- 
cause they are 'A' effects, i.e., the fractional change in 
energy is inversely proportional to the number of elec- 
trons in the system. The precision of the calculation 
must therefore be sufficient to resolve this energy change 
amid the statistical noise. The system must also be large 
enough to give a good description of the infinite solid, 
which is a severe constraint for small band gap materials 
such as our test material, silicon. So far only a few DMC 
calculations of excitation energies in solids have been re- 
ported. In particular, we note the 8-atom simulation 
cell calculations of an energy gap in a molecular nitrogen 
solidtil and the T^y — > X\ c and Yi u — > X\ c excitations in 
carbon diamondu. The excitations in these calculations 



were indirect in reciprocal space, so that the excited state 
is orthogonal to the ground state hy translational sym- 
metry. In previous work on solidsoO it has been assumed 
that the standard DMC method would give good results 
only for the lowest energy state of each symmetry. How- 
ever, in this paper we will show that the DMC method 
can be applied successfully to a wide range of excitations 
in solids. 

We have chosen silicon for our study because (i) DMC 
gives a good account of the ground stateQ, (ii) electron 
correlations significantly affect the band energies, and 
(iii) results from other [Calculational methods, such as 
the LDA, HF0, and GVFtrO approximations are well es- 
tablished, and a large amount of experimental data also 
exists. In this paper we explore the limits of the DMC 
method by calculating excitations which are direct and 
indirect in reciprocal space, and including several excita- 
tions at each wave vector. 

It is important to distinguish between different types of 
excitation energy. In quasiparticle theory the quasipar- 
ticles correspond to the poles of the one-particle Green 
function and are equal to the energies for adding an elec- 
tron to the system or subtracting one from it. The quasi- 
particle energies have both real and imaginary parts, the 
latter giving the quasiparticle lifetime. These energies 
are measured in photoemission and inverse photoemis- 
sion experiments. For the minimum gap the imaginary 
part of the quasiparticle energy is zero and the quasipar- 
ticle has an infinite lifetime. In this case the quasiparticle 
energy gap can be written as E g = Ej^ + i +£jv_i — 2En, 
where -Ejv+i, -Ejv-i, and En are the ground state total 
energies of the N+l, N — l and N electron systems. This 
energy gap is accessible within DMC methods, but we do 
not consider it here. In an optical absorption experiment 
a different process occurs in which an electron is excited 
from the valence to the conduction band. In this case 
an exciton is formed and the lowest excitation energy is 
smaller than E g by the exciton binding energy. In the 
calculations reported here we create excitonic states by 
exciting electrons from a valence band state into a con- 
duction band state. Although the exciton binding energy 
is artificially increased by the finite size of our simula- 
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tion cell, it is still small (about 0.1 eV) and therefore 
our results are comparable with theoretical quasiparticle 
energies and experimental photoemission data. We re- 
mark that excitonic properties are of significant interest 
in their own right and the ability to calculate both elec- 
tron addition/subtraction energies and excitonic energies 
is a significant advantage of the DMC method. 

This paper is organised as follows. In Sec. II we pro- 
vide a brief description of our calculational method. Sec- 
tion III is a detailed discussion of the results of our study 
demonstrating the viability of the DMC method for ac- 
curately determining the low lying excitation energies in 
solids. We conclude with a summary in Section IV. 



II. CALCULATIONAL METHODS 

In the DMC methodS! imaginary time evolution of 
the Schrodinger equation is used to evolve an ensemble 
of 3./V-dimensional electronic configurations towards the 
ground state. Importance sampling is incorporated via 
a guiding wave function, $. To make the calculations 
tractable we use the fixed node approximation, in which 
the nodal surface of the wave function is constrained to 
equal that of $. The fixed- node DMC method gener- 
ates the distribution Q*i>, where if? is the best (lowest 
energy) wave function with the same nodes as The 
accuracy of the fixed node approximation can be tested 
on small systems and normally leads to very satisfactory 
results! We also use the short-time approximation for 
the Green's function, whose effect can be tested and made 
very small. We used a time step-pf 0.015 a.u., which gives 
small time-step errors in siliconEl. The average number of 
configurations in the ensemble was 384, and between 1220 
and 2125 moves of all the electrons in all the configura- 
tions were attempted, except for the ground state where 
3848 moves were attempted. In all cases the acceptance- 
rejection ratio was greater than 99.7%. 

We used a fee simulation cell containing 16 silicon 
atoms, employing periodic boundary conditions to reduce 
the finite size effects. The Si 4+ ions were represented by 
a norm-conserving non-local LDA pseudopotential, and 
the non-local epergy was evaluated using the "locality 
approximation"!!!]. The non-local potential was sampled 
using the techniques of Fahy et al.H3 Our guiding wave 
functions are of the Slater- Jastrow type: 
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where there are N electrons in the simulation cell, \ is 
a one-body function, u is a two-body correlation factor 
which depends on the relative spins of the two electrons, 
and and are Slater determinants of up- and down- 
spin single-particle orbitals. We used a Fourier series 
expansion for the \ function which was constrained to 



have the full symmetry of the diamond structure and con- 
tains 6 free parameters. The parallel- and antiparallel- 
spin u functions were constrained to obey the cusp con- 
siditionsEJ and contained a polynomial part with 11 free 
parameters, as described in Ref. [l5|. The guiding wave 
function contained a total of 28 parameters, whose op- 
timal values were obtained by minimizing the variance 
of the energy using- ID 5 statistically independent elec- 
tron configurationaLS'Ej, which were regenerated several 
times during the minimization procedure. The optimal 
parameter values were obtained for the ground state wave 
function and were used for the excited states as well. Be- 
cause the parameters occur only in the nodeless Jastrow 
factor this procedure does not bias the DMC excitation 
energies, which depend only on the nodal surfaces of the 
guiding wave functions. 

Our calculations are for many-body Bloch wave func- 
tions having definite values of the crystal momentum or 
wave vector, k. The Slater determinants were formed 
from the LDA orbitals calculated at the T-point of the 
Brillouin zone of the simulation cell, which unfolds_to the 
r, X and L points of the primitive Brillouin zoneEJ. The 
wave vector of a determinant of these orbitals is equal to 
one of the r, X or L wave vectors. The determinant for 
the ground state guiding wave function was constructed 
from the valence band orbitals at the T, X and L points, 
and has k=0. The excited state guiding wave functions 
were formed by replacing an orbital in either the up- or 
down-spin determinants of the ground state wave func- 
tion by a conduction band orbital. 

The computational demands of DMC calculations are 
such that presently it is not feasible to calculate excita- 
tion energies with cells larger than our 16 atom cell, and 
therefore it is important to consider the finite size effects. 
The finite size effects can be divided into "independent- 
particle finite size effects" (IPFSE), which can be mod- 
elled by LDA calculations, and "Coulomb finite size 
effects" (CFSE), which arise from the explicit use of 
Coulomb interactions between the particles. The HF, 
LDA and DMC ground state energies for the 16 atom 
cell are -103.67, -106.61 and -107.31(1) eV per atom, 
respectively. We estimate finite size errors in the ground 
state calculations by performing calculations on simu- 
lation cells containing 250 atoms (HF and LDA) and 
128 atoms (DMC), giving corrections of -0.48, -1.31 
and —1.02(3) eV per atom, respectively. The IPFSE are 
therefore quite large, being about 21 eV for the 16 atom 
cell used here. However, LDA calculations performed 
for a number of simulation cell sizes show that the can- 
cellation between the errors in the ground and excited 
states is so good that the resulting finite size errors in 
the LDA excitation energies of the 16 atom cell are less 
than 0.1 eV. We expect a similar cancellation of errors in 
the DMC cancellations. In our DMC simulations we have 
also used a new formulation of the electron-electron inter- 
action within periodic boundary conditionals which gives 
much smaller CFSE than the standard Ewald forrrca. 
However, the CFSE also tend to cancel in the excitation 
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energies and we obtained almost identical results using 
our new interaction and the Ewald interaction. 

Each of the guiding wave functions we use contains 
some spin contamination, i.e., they are not eigenstates of 
S 2 but are admixtures of different spin components. We 
have investigated the effect of this spin contamination by 
calculating the T25' —>Xi c excitation energy using a spin 
contaminated single determinant and a two-determinant 
spin-singlet wave functioned. These calculations gave en- 
ergies differing by less than the statistical noise of 0.2 
eV. The exciton binding energy in our calculations is en- 
hanced because the exciton is confined to the simulation 
cell. Following Ref. || we use the Mott-Wannier formula 
to approximate the binding energy of the localized exci- 
ton, giving 0.1 eV, which is smaller than the statistical 
noise in our calculations. Finally we investigated the ef- 
fect on the DMC energy of using single-particle orbitals 
which had been relaxed by performing LDA calculations 
in the presence of the excitation. The resulting changes 
in the excitation energies were smaller than the statistical 
noise. 



III. RESULTS AND ANALYSIS 

In Table I we give our results for the 27 excitations 
studied, together with HF, GW, and LDA data. The 
characteristics of the HF, GW, and LDA excitation en- 
ergies are well known. The GW approximation gives ex- 
tremely good excitation energies for weakly correlated 
systems such as silicon, while the LDA excitation ener- 
gies are too small by 0.7-1.0 eV, and the HF excitation 
energies are much too large. The agreement between the 
DMC and GW excitation energies is good for the low 
energy excitations, but poorer for the higher energy ex- 
citations. The percentage, ay, of the correlation energy 
retrieved by our DMC calculation for the state formed 
by exciting an electron from single-particle orbital i to j 
is 

aij = t ^m + at) ^ c x ioo . (2) 

K lJ — E lJ 4- E° 

^ exact ^HF ' 11 c 

The values of E^) MC and E^ F are given in Table I. We 
estimate the correlation energy for the ground state, E®, 
to be about -60 eV per simulation cell, and we estimate 
the fraction of the ground state correlation energy re- 
trieved by the DMC calculation, a , to be in the range 
90-100%. For the purposes of this comparison we use 
the GW energies from Ref. [ll] for the E l J xact because of 
the incompleteness and uncertainty of the experimental 
data. |-We note that the various GW calculations for sili- 
conEliiil are in excellent agreement with one another and 
are also in good agreement with the available experimen- 
tal data. The experimental data are reported as band en- 
ergies, so that to form the excitation energies we have to 
take differences between the experimental values, which 
adds to the uncertainties. When we present our results 



for band energies we will compare with the experimental 
data, which is given in Table II. 

The values of ay given by this analysis slowly de- 
crease with increasing excitation energy, so that for the 
largest excitation energies the a y - are 2-3 percentage 
points smaller than for the smallest excitation energies. 
This conclusion is insensitive to the values of E® and ag. 
This suggests the following rationale for our results. The 
HF excitation energies are much too large because the 
correlation energy is neglected. The fraction of the cor- 
relation energy retrieved by our DMC calculations slowly 
decreases with increasing excitation energy. However, be- 
cause the contribution of the correlation increases rapidly 
with increasing excitation energy the DMC excitation en- 
ergies are somewhat too large. This analysis indicates 
that the residual errors in the DMC excitation energies 
are mostly due to the errors in the nodal surfaces of the 
excited state guiding wave functions rather than in the 
ground state, and that the quality of the nodal surfaces 
falls with increasing excitation energy. 

To study whether the DMC method works for direct 
excitations we considered pairs of excitations where a 
single-particle orbital with a particular wave vector is re- 
moved from the determinant and replaced (i) by a higher 
energy orbital at the same wave vector and (ii) by an 
equivalent higher energy orbital at a different wave vec- 
tor. For example, when calculating the A4 — » X\ c exci- 
tation energy, we replace an X4 orbital at a particular 
X-point by an X\ c orbital at the same X-point to give 
a direct excitation, while for the X4 — >X^ C excitation we 
replace the X4 orbital with an X\ c orbital from a differ- 
ent X-point, giving an indirect excitation. We have in- 
vestigated three such direct-indirect excitation pairs and 
found reasonable consistency between the results. We 
have calculated a total of 6 direct (r-point) excitations 
and the level of accuracy is indistinguishable from that 
for indirect excitations, showing that the DMC method 
can be applied to direct excitations. In addition we have 
successfully calculated a total of 10 excitation energies at 
the X-point and 11 at the L-point, which demonstrates 
that the DMC method works for higher excitations as 
well. 

The 27 entries in Table I correspond to transitions 
between 7 valence and 5 conduction band energy lev- 
els. To obtain DMC band energies, e,-, we performed 
a least squares fit to the DMC data by minimizing 
T,[Ep MC — (ej — tj)] 2 with respect to the e,, where the 
sum is over the 27 excitation energies, Ep MC , listed in 
Table I. The resulting band energies are given in Table II 
together with other theoretical data and with experimen- 
tal data. For greater clarity the DMC band energies are 
plotted in Fig. 1 together with an empirical pseudopo- 
tential band structureEj which is in good agreement with 
the available experimental data. The energies at the top 
of the valence band have been aligned. The DMC band 
energies are very much better than the HF values because 
of the inclusion of correlation effects. For the lower part 
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of the valence band the DMC energies lie consistently 
1 to 1.5 eV below the empirical pseudopotential, GW, 
and experimental data. The DMC band energies around 
the gap region are in good agreement with the empirical 
pseudopotential, GW, and experimental data. 

The success of our excitation energy calculations is 
very encouraging. To appreciate why our calculations 
are so successful we must consider the fixed-node DMC 
algorithm in more detail. First we consider ground state 
calculations. The fixed- node DMC energy is a variational 
bound on the exact ground state energy, and if the nodal 
surface of the guiding wave function is exact then the 
resulting fixed-node DMC energy is also exactEl The 
exact ground state wave function "tiles the configuration 
space" , which means thai, all nodal pockets are related by 
permutation symmetrycj. Normally one chooses guiding 
wave functions which also have the tiling property and 
consequently the DMC simulation may be performed in 
any subset of the nodal pockets. 

For excited states the situation is somewhat different. 
If the fixed-node constraint could be removed the algo- 
rithm would, in the limit of a long simulation, give the 
ground state energy. However, the key point is that the 
fixed-node approximation prevents this collapse and al- 
lows us to obtain good approximations for excited state 
energies. One can readily show that if the nodal surface of 
the guiding wave function is exactly that of the nth eigen- 
state, then the fixed-node DMC procedure gives the exact 
energy of the nth eigenstate, regardless of whether the 
guiding wave function overlaps with lower energy states. 
However, if we use a guiding wave function that does not 
have the exact nodal surface of the excited state we are 
modeling, the resulting DMC energy may not be above 
the exact energy of that state because of the possibility 
of mixing in lower energy states with the same symmetry. 

Another complication arises from the fact that the 
nodal pockets of approximate excited state wave func- 
tions may be inequivalent. Therefore, as the DMC simu- 
lation proceeds, the population of configurations in some 
nodal pockets will dominate over others, and in principle 
the results could depend on which of the pockets were 
occupied at the start of the simulation. An illustration 
of this type of behaviour for the first excited state of an 
electron in an infinite square well is given on page 186 of 
Ref. ||, which shows the DMC energy decreasing linearly 
with the magnitude of the error in the nodal position. 

Although there are significant differences between 
ground and excited state fixed-node DMC calculations, 
the criterion for obtaining good energies is the same, i.e., 
the nodal surface of the guiding wave function must be 
of good quality. Our results demonstrate that the nodal 
surfaces from determinants of LDA orbitals are fully ad- 
equate for calculating excitation energies around the gap 
region in silicon, but give poorer results for higher ex- 
citation energies. It will therefore be necessary to opti- 
mize the nodal surfaces of the guiding wave functions to 
obtain more accurate estimates of the higher excitation 
energies. We note that a DMC calculation by Grimes et 



al. for an excited state of the hydrogen molecule with 
the same symmetry as the ground state also gave a good 
excitation energy. This result lends further support to 
our contention that the fixed-node DMC method can be 
applied to a wide range of excited states. 

IV. SUMMARY 

In conclusion, the DMC method is a stable and accu- 
rate method for calculating low lying excitation energies 
in solids. The fixed- node approximation works to our 
advantage by preventing collapse to lower energy states. 
The accuracy of the excited state energies is determined 
by the quality of the nodal surfaces of the guiding wave 
functions. We have obtained good values for the low lying 
excitation energies, including direct and indirect excita- 
tions, and including several excitations at each wave vec- 
tor, which indicates that the nodal surfaces of our guid- 
ing wave functions are of good quality. We have demon- 
strated that DMC calculations can be used to calculate 
direct excitations, which is important because it allows 
us to obtain excitation energies when there is no un- 
derlying translational symmetry. The fixed-node DMC 
method provides a unified framework for calculating ac- 
curate ground and excited state energies. 
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TABLE I. Excitation energies in eV calculated with the 
DMC, HF, GW, and LDA methods. The wave vector of the 
excited state wave function is denoted by k. The statistical 
error bars on the DMC energies are ±0.2 eV. 
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TABLE II. Bp«d energies of silicon in eV. a- This work, 
b- ReflU c- Ref.E-!J. d- Ref.1 2 -!!. e- From the compilation given 
in Ref.E-1 The statistical error bars on the DMC energies are 
±0.2 eV. 



Excitation 


k 


DMC 


HF 


- 


GW 


Ll 


LDA 


Tiu 




r 


18.05 


27.9 


15.84 


15.14 


Fiv 


— >ri5 


r 


17.38 


26.9 


15.31 


14.50 


Tit, 


— * Lie 


L 


16.27 


25.4 


14.14 


13.39 


Tit, 


— > Yic 


X 


14.93 


24.2 


13.38 


12.58 


r 25 < 


->L S 


L 


4.75 


8. 


7 


4.25 


3.31 


r 25 < 


^Tis 


r 


3.82 


8. 





3.36 


2.55 


r 25 ' 


— * Lie 


L 


2.35 


6.5 


2.19 


1.44 


r 25 < 


— > Yic 


X 


1.34 


5.3 


1.43 


0.63 


X\ v 


^r 15 


X 


12.42 


20.5 


11.31 


10.36 


X\ v 




L 


10.37 


17. 


8 


9.38 


8.44 


X 4 - 


^Fl5 


X 


6.91 


12.7 


6.29 


5.39 


Y 4 - 




X 


7.92 


13.7 


6.82 


6.03 


Y 4 - 


-> Lie 


L 


5.71 


11.2 


5.12 


4.28 


Y 4 - 


-> Yic 


r 


5.12 


10.0 


4.36 


3.47 


Y 4 - 


~* -Ale 


L 


4.86 


10.0 


4.36 


3.47 


L 2 , 


^L% 


X 


15.60 


24.1 


13.95 


12.92 


L 2 , 


^ris 


L 


14.75 


23.4 


13.06 


12.16 


Li v 




X 


12.29 


19. 


S 


11.39 


10.30 



5 




FIG. 1. The DMC band structure (filled circles with error 
bars) . As a guide to the eye we also show empirical pseudopo- 
tential dataLJ (solid lines). 
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